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ABSTRACT 


An overview is given of multivariate, wide-sense, sta- 
tionary, stochastic process, linear prediction theory with 
emphasis on the mean square error of prediction based on 
a finite past. Examples of several different types of proc- 
esses are examined with specific formulas for the pre- 
diction errors obtained in some cases. A discussion is 
given concerning the effect on the prediction error of basing 
the prediction on a subprocess of the given process instead 
of on the entire original process. 
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LINEAR PREDICTION ON A FINITE PAST 
OF A MULTIVARIATE STATIONARY PROCESS 

by 

Ray G. Langebartel* 

Goddard Space Flight Center 


INTRODUCTION 

Linear prediction theory for one-dimensional stationary stochastic processes ("time series") 
has been studied extensively (for the essential aspects see, e.g., Reference 1). Less attention has 
been paid to multivariate processes but a body of theory is now at hand (References 2 to 5). Pre- 
diction on a finite past in a multivariate process with discrete parameters is not dwelt on in these 
works, however. In the present paper this aspect of the theory is considered with special emphasis 
on the mean square error of prediction. In certain quite special multivariate processes fairly 
explicit results are obtained, but the problem for more general processes is difficult. 


BASIC CONCEPTS 

Let X m , where m = 0, ±1, and ±2, * * * , represent a multivariate discrete stochastic process of 
dimension M : 


a column matrix of M rows. The x m k may be complex numbers, in which case we designate the con- 
jugate of x m k by x m k and the conjugate of the matrix X m by X^, the matrix of the conjugates. The 
transpose of a matrix will be designated by T in this manner: X m T . The expected value of a multi- 
variate process is taken as the matrix of the expected values of the elements: 

E k}' 

E k 2 } • 

E k H }_ 

*GSFC consultant from the University of Illinois. 
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The covariance matrix (or correlation matrix) E{x m X n T } is defined accordingly: 
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In the event that each of the elements of this matrix, i.e., the various covariances of the x m k proc- 
esses, depends on m and n only through their difference the multivariate process is called sta- 
tionary in the wide sense (or weakly stationary) and we write 


E k X „ T } S Vn) " R m-n • 

It is a straightforward computation to show that the covariance matrix of a wide-sense stationary 
multivariate process satisfies the relation 



R (-m) 


( 1 ) 


It is often the case in practice that the elements of the ftiatrices X n and R n are of different 
magnitudes and it is desirable to scale them to some norm. The elements of X n can be scaled by 
premultiplying X n by a diagonal constant matrix A, 


to give 



0 0 • • • o 

k 2 0 ••• o 

0 k, ... o 
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k 2 X n 2 


k M X " 
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For the scaled covariance matrix R m we then have 


K s E{X m+n X n T } = E{AX m+n (AX n ) T } = E{AX m+n X n T A T } = AE{x m+n X n T } A T = AR m A. 

It is customary to choose the elements of A so that the elements of R 0 down the main diagonal are 
all equal to unity. In the case of real matrices this means that A is 



where R 0 lj = E{xj xj } , the elements of R 0 . For M - 2 the scaled covariance matrix is, when written 
out, 
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One way to investigate the nature of aprocessis to treat its Fourier representation. The wide- 
sense stationary process X m has the spectral distribution function (square matrix) F(X) represent- 
ing the covariance matrix, 


V, = 1 


1/2 

1/2 


R,_ x = | e 2TTi*.n dF (\) 


with a corresponding representation of the process itself, 


r 1/2 
■* — 1 / 2 


X_ = | e 27 ^ dY(A-) 

- 1/2 


( 2 ) 


(3) 
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where Y(X) is a column matrix (Reference 1, p. 596; Reference 6, p. 45). There must, of course, 
be some relation between F(X) and Y(X) for these representations to hold. In terms of the spectral 
density function F' (X) , this relation is 


E{dY(X)dY(>) T } " F’ (X) S(X - /x) dX d/x 


(4) 


where S(X) is the Dirac delta function. The Y(X) process is spoken of as having "orthogonal in- 
crements." The reason for the form of Equation 4 is indicated by the following heuristic calculation 


f r 1/2 r 1/2 

R (m) = E{x m+n x7} - e 2iri(m+n)X dY(X) dY( M ) T 

[J-l/2 J-l/2 
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which is independent of n and thereby confirms wide-sense stationarity. In a particular example 
F(\) may well have jumps so that F' (\) would not exist in the classical sense (although it would in 
the Schwartz distribution theory sense) and this is the reason that in much of the literature these 
Fourier representations are presented in the Stieltjes form. 

If the process is real the representations of Equations 2 and 3 can be put in a real form. Let 
dY = dYj + idY 2 , SO that 

/* 1/ 2 

x m - ^(cos 2mnA.) dYj - (sin 277mX) dY 2 J + i [(cos ^nA) dY 2 + (sin 2rrmk) dYj ] . 

J-l/2 

This is evidently real if dY 2 is an odd function of k and dY l is even (i.e., if Y 2 ' (k) is even and Y^A.) 
is odd), in which case the real part of the integrand is an even function of k and only the half 


4 



interval (0, 1/2) is required: 


2 



[/cos 2 mnk) dYj 


(sin 2rrm \ ) dY 2 ] . 


Let 


2dYj = dU 


and 


- 2dY 2 ^ dV 

and, by substitution, obtain the representation for real values of X m 

r 1/2 

X m - I (cos 2miA) dU(X) + ( s in 2nnA) dV(X) . ( 5 ) 

Jo 

A corresponding treatment on Equations 2 and 4 leads to the companion formulas 

r 1/2 

R (m) = (cos 27rX.m) dG(\) , (6) 

E{dU(\) dU ( T M) } = E{dV(\) dV^} = G' (\)8(\- M )d\d M 

E{dU(\)dV ( T M) } = E{dV(X) dU ( T M) } = 0 . 

It is natural to introduce orthogonality terminology in view of the relations satisfied by dY, du, 
and dv. In particular, we say that X and Y (column matrices with random variables as elements) 
are orthogonal if 



e{xY t } = 0 . (8) 

Furthermore, X and Y are orthonormal if in addition to Equation 8 they satisfy 

e{»?} = e{yyt} = I (9) 

N 

where I is the identity matrix. An expression of the form £ A n X n where X 1? X 2 , • • ■ X N are 

n — 1 

multivariates of dimension M and , ■ ■ • A N are MxM matrices of constants is called a linear 
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combination of the X n . The collection of all linear combinations of x l9 * • ■ X N forms the linear 
manifold generated by X 19 * * • X N . The orthogonal projection of Y on a linear manifold is Y p if Y p 
belongs to the manifold and Y - Y p is orthogonal to every element of the manifold. 


LINEAR PREDICTION 


Suppose in a wide-sense, stationary, multivariate, stochastic process X m we wish to predict 
X N+1 given the values of x 0 , X 1? X 2 , * * • X N , and that this predicted X N+1 can be represented as a 
linear combination of x 0 , • • • X N . One method of attack is to introduce H N+1 , which is to be the unit 
normal of the manifold spanned by X 0 , * * * X N , and to represent it as a linear combination of 

^0 ? 9 ^N + 1 ‘ 


Then 


and 


Hn +1 = 22 B r X r • 

r= 0 

{ N + l 

x s £vv 

{ N + l ~j 

22 

N+l 

22 e{x s x7} §7 


N + l 




R B t 

S- r r 


for 0 < S < N 


r v n+ i 

I - e{^ + l ~ N + 1 } “ E | E r x r a N+1 


n Ti 

22 B r E K^ tl } 

r= 0 

En+i eK+i =n T +i } 

N+l 

®N+1 y ^ ^N + l-r ' 
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These two results can be combined to give the matrix equation for C r = b7b n+1 , provided B N+1 is 
nonsingular. 


N+l 

L 

r-n 


*s-r c r 


s 


s,N+l 


I , 


for 0 = s ^ N + l , 


(10) 


where S , is the Kronecker delta. But 


B N+1 ^N+l “ / Vl X r B N+1 ®N+1 + 


r= 0 


so if we use ¥ as the predicted value for x N+1 then the mean square error is 


E {(V, - *) ( x n +1 - *) T } ^ e{b n ;\ h n+1 ( B - + \ H^)*} 

E K+1 Hn+! ®n+i b n +7 } 


Vl E{eN +1 V-i } B i 


IT 
'N + l 


13 ”1 o ~1T 
^N+l °N+1 


C 

^N+l * 


The best predicted value in the mean-square sense for x N+1 is ¥. This may be seen intuitively 
on geometric grounds by observing that ¥ belongs to the manifold spanned by X 0 , • • • x N and 
X N+ i -¥ = ^n+i is Perpendicular to the manifold, so that¥ is the orthogonal projection of X N+1 
on the manifold. Thus, ¥ is the closest vector in the manifold to X N+1 . 

This line of reasoning has its counterpart in the following development. Let 0 be any vector 
in the manifold. Then 


E {( X n + i -®)( X n + i -®") T } = E {[( x n + i -'P) + ( , P -®)][(X n+1 -*) + (»■ -0)]} 

= E {( X rm -’•’)( x nh-' I ') T } +E {( X n« ~ }+ - 8) } + E {(T~ 

= E {( X n + i -«)( x N+ i r ^} + E {(V-0)(>^ r 0 ) T }- 
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If these various nonzero expectation matrices are nonsingular then they are Hermitian symmetric 
and positive definite, so that for an arbitrary vector z 


zTe {( X n«-®)( X n + i-®) T } Z = zT E {( X n + i - 'P ) ( X N + 1 - ^) T } Z + Z t e{('P-0)(W-0) t }z 

(where each term is real and positive), and this has its minimum at 0 = w. 

Consequently, the multivariate linear least- square prediction of the x N+1 problem given 
X 0 , • * • X N becomes the problem of solving matrix Equation 10 for c r where the R s __ r are regarded 
as known. The predicted value is 


N 

* = - ^B N V lBr X r 

r=0 


and the mean square error is - B N+ 1 1 . It should be remembered that the X m process is 
assumed stationary in the wide sense. 

If the prediction were to be made on the basis of the entire (infinite) past, that is, on the basis 
of X m where -°°<m< N, then Equation 10 becomes in effect an infinite number of equations with an 
infinite number of unknowns and as such is unruly. In the one-dimensional case it has been found 
advantageous to employ the Fourier representations of Equations 2 and 3, in which case the mean 
square error turns out to be 


e 


1 r 

2 I log F' (X) dX 
*'- 1/2 


(ii) 


where F' (X) is the (one-dimensional) spectral density function of the process (Reference 1, p. 577). 
Rosenblatt (Reference 2), among others, extended the theory of prediction on an infinite past to the 
multivariate case but he discovered that the prediction error does not turn out to have the simple 
form of Equation 11 (where F' (X) is the matrix multivariate spectral density function) except in 
rather special cases. The difficulty arises from the noncommutativity of matrices in general. 

The above results for prediction based on a finite past can also be obtained through use of the 
Fourier representation of Reference 2. Suppose that the spectral distribution function F(X) (an 
MxM matrix) is absolutely continuous so that there is a spectral density matrix f' (X) which we 
assume to be continuous and nonsingular for all values of X. Let the spectral density matrix be taken 
as 


F' (X) 




( 12 ) 
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where 


' N+l 

L 

\ k~0 


B T Z N + l-k 


is nonsingular for |z| ^ 1. Thus the "Hermitian square” of the matrix H (k) is F' (A.): 

HH^ - F' (k) 


where 


H(M 



(13) 


Introduce this "square root" matrix into the Fourier representation of the X m process in Equation 3 
in the fashion 


r 1/2 

x m = e 277inA H(A)d <J>(\) (14) 

1 / 2 

where d$ is a column matrix satisfying 


E{dO(X)d<I>(/x) T } = S(k-fJL)Idkd/jL ; 


then 


E{X m+n Xj} 


r 1/2 

r 1/2 

J -1/2 

" ~l/2 

r 1/2 

e 27rin\ 

■'-1/2 



e 27ri(m+n)A.-27,im K g(A. - p.) I d\ d,aH(,u.) T 


H(X)H(A.) T dA. 


• 1/2 

e 277inX p/ , 

- 1/2 


which agrees with Equation 2. The justification for these manipulations (including the existence 
of d<&(\)) is discussed by Rosenblatt (Reference 2). Now we make use of the special form we have 
assumed for H (A). Define H n by 


^n+N+l 


N+l 



k = 0 
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then 


-n+N+l 


N+i r 1/2 

^ B k J e 27ri ( n+k A H(X) d<I>(X) 

k=0 ~ 1/2 

/•1/2 N+l / N+l 

J e 277in\ B k e 277 **- ( ^ B k 


d<D(\) 


r 1/2 

= J ' 

*'- 1/2 


e 27TinX 


which shows that the d<t> process represents E n+N+1 just asHdfl) does X n (see Equations 3 and 14). 
The above apparatus was set up precisely so that E n would be orthogonal to X m for n > m and so that 
{ E n }would be an orthonormal sequence. This gives the same orthogonal projection and linear 
manifold situation that was arrived at in the earlier treatment except for the greater generality 
now in that we are dealing with the sequence X n , X n+1 , • • • X n+N+1 for any n instead of the one se- 
quence X 0 , Xj , • • • X N+1 . To show the orthogonality of E n with X m , 


r r 1/2 /M/2 

E{HnXj} = E< e 2 7 /i( n -N-l)\ d<I>(X) j e" 277 ^ (H(/Li) d4><») 

l/- 1/2 J-l/2 


[ 1/2 fl/2 

e 277i(n-N-lA-277i m/ x dd>(» T H(/T)T" 

y- 1/2 ^ — 1/2 


"1/2 •'-1/2 


fl/2 fl/2 

= e 2^i(n-N-l)X-27/im M /J.) I dXd/U-HCM) 1 

J~1/ J •/-I/O 


-1/2 ^ - 1/2 


r 1/2 

= e 27/i(n-„,-N-i)\^yr d\ 

•'“l/ 2 


ri/2 

= J e 

" “ 1/2 


1/2 / N+l 

e 27ri(n-m-N- 1)A. | \ * 


B T g-2/rikX ^ 


“ 1 


2^1 j* Zn ‘"’" 1 ^ B k T Z N+1 - k j dZ 

0 for n - m = 1 , 

(B7 + i) 1 for n = m 
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The final integral, which is a contour integral around the unit circle, was obtained by making the 
transformation e 277 ^ = z. The original defining assumption on H(\) in effect said that 


N + l 

L 


k=0 


B k T Z N + l-k 


is analytic on and inside the unit circle. Consequently, ifn-m-l^o, the entire integrand is analytic 
on and within the contour so that the integral vanishes. On the other hand, if n = m then the factor 
z n- m -i Yi&s a simple pole at the origin and an immediate residue calculation gives the integral’s 
value as (b^) -1 . 

To show the orthonormality of E n : 


f r 1/2 i * 1/2 

E{s n sj} = E*s d$( M ) T 

ly-1/2 J-l/2 

r 1/2 pi/ 2 

J -I/O J- 


5 27ri(n-N-l)\-277i(m-N-l) M - fi) I dk dfJL 


-1/2 ^“1/2 


r 1/2 

e 277i(n-m)\ j dk 

•'- 1/2 


The matrix Equation 10 is also easily derived by the Fourier method. 



N+l 

L 

k = n 


B k R 


T 

n+s~k 


p 1/2 / N + l 

J e -277i(s+n)X / ^ §T e -277 ik\j d\ 

“1/2 \ v-n / 


N+l 

£ B k R k-s-n 

k = 0 


(by Equation 1) 


N+l r 1/2 

B k e 2iri(k-s-n)\ F' (\) d\ 

W J - 1/2 

2^1 zN ' s '" zN+1 ' k ^ dz 


B T ” 1 
d n+i 


for N - s - n ^ 0 , 
for N-s-n = -1 . 
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Upon taking the conjugate transpose we obtain 


N+l 

E 

k=0 


■p o T — D-l S 

^n+s-k D k °N+1 s.N + l-n 


for s ^ N + 1 - n , 


which agrees with Equation 10 when n = 0. 


PREDICTION ERROR MATRIX EQUATION 

The determination of the mean square error matrix for linear least-square prediction of x N+1 , 
given the values of x 0 , • • • where X m is a wide-sense stationary process, is equivalent to solving 
matrix Equation 10 for C~* l9 which is the mean square error. This equation, while appearing to be 
a fairly general linear equation, is actually a bit special in that the coefficient matrices R k are 
subject to the relation of Equation 1. In itself this does not appear to be of much help when it comes 
to a straightforward computation of c N “* from the simultaneous linear set, but further conditions 
on R k may aid materially. In particular, we shall show that if R k is Hermitian symmetric then the 
simultaneous set breaks up into two sets, each of approximately half the order of the original set. 

Let k = n + l “ s in Equation 10. Then the equation is 


N+l 

2^ Wk-r C r = S k0 I for 0 = k I N + 1 . ( 15) 

r= 0 

Written out and upon application of Equation 1 this becomes 

Vl C 0 + Vl + Vl C 2 + - +R 2^1 +R l^ +R 0^1 * 1 

^ c 0 + Vi c i + V 2 c 2 + ••• +R 1 ^- 1 +R 0 qv + R ^ +1 = 0 
Vi c o + + Va c 2 + •" + R o<Vi + R 7 ^ +r 7<^ +1 - 0 


R 2 C 0 + Rj Cj + Rfl C 2 + • • • + R,7_3 + Rp^L 2 + ^n-i *-n+i ® 

RjCo + RoC, + ^c 2 + - +^T 1 q,+^? r q l+1 = o 

R 0 c 0 +5/C, +r/c 2 + + ^q, +1 = o . 
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If we let R represent the matrix of coefficient matrices, C the matrix of the matrix unknowns, and 
J the column matrix of matrices on the right, then this system can be written as 

RC = J . 

Suppose the matrices R m are all Hermitian symmetric. Then the matrix R takes on a good deal of 
symmetry, both centrosymmetry (symmetry about the center point, i.e., the point of intersection 
of the two main diagonals) and persymmetry (all elements of any one line perpendicular to the 
major diagonal alike). We shall treat only the case where N is even; the case for odd values of N 
can be handled similarly. Introduce the matrix D of the identity and zero submatrices: 


and 


“I 

0 

0 

— 

0 

0 

r 

0 

I 

0 

— 

0 

I 

0 

0 

9 

I. 

•II* 

I 

9 

9 

0 

0 

O' 

:??:■ 

■ I 

0 

0 

0 

0 

0 

— 

0 

i 

0 

_0 

0 

0 


0 

0 

i 


I 0 0 0 0 -I 

0 I 0 0 -I 0 

0 0 I. .t-I 0 0 

: : : 'I-I : : : 

: : • *0 I* : : : 

0 0 0 * T 0 0 

0 0 0 0 I 0 

0 0 0 0 0 I 


The inverse matrix D~ 1 is to be understood in the sense that DD 1 is to be the matrix with the 
identity matrices I down the main diagonal and zero matrices elsewhere. 

Now 


R N + 1 

Kn 

r n-i 

R. 

£N + 1 

R, - R, 

5N ^N+l 

R, - R, 

jN~l 5N + 2 

R 2 ** R N - 1 

R 1 ' r n 

R 0 R N +1 

R n 

r n-i 

R N -2 

Rj 

2 N 

R 1 " R 1 

R, - R, 

jN "2 ^N+l 

R 1 “ R N “2 

R 0 " R N -1 

r i “ r n 

r n~i 

R N ”2 

R N -3 

Ri 

2 n “ 1 

r i “ R i 

jN-2 ^N-l 

R. - R, 

jN “3 5N 

R 0 ~ R N -3 

R 1 " R N “2 

R 2 " R N -1 

R, 

^N + l 

R 1 

R 1 

R 1 

R 0 ~ R 1 

Rj - R 2 

R, -R t 

^N“2 ^N-l 

R 1 _R 1 
An-1 An 

R 1 “ R 1 
2 N 2 w+l 

R, + R, 

jN+l 

R 1 +R 1 
jN“l 

Rj +R. 

jN~2 ^N-I 

R o +R i 

0 

0 

0 

0 

0 

R, + R, 

£n -1 jN + 2 

R, + R, 

^N“2 ^N+l 

R 1 +R l 

^“3 Jn 

R i + R 2 

0 

0 

0 

0 

0 

R 2 + R,^ 

R 1 ^N-2 

R 0 + R N -3 

r i + *Ri 

^N-2 ^N~l 

0 

0 

0 

0 

0 

R 1 +r n 

R 0 + R N -1 

R 1 + R N -2 

• • ■ R, + R, 

jN“l jN 

0 

0 

0 

0 

0 

R 0 + R N + 1 

R 1 +r n 

R 2 +R N -1 

R, +R t 

■jN ^N + l 

0 

0 

0 

0 

0 

_ 
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DC 


+ ^+1 
C i 

c 2 + q*-! 

C 1 + C 1 
5 n £n + 1 

C 1 
iN + i 

$ 1-1 
^ + 1 


and 


D T J 


I 

0 

9 

o 

9 

o 

0 

1 


The equation RC - J is equivalent to D T RD 1 DC = D T J and in view of the three matrices just cal- 
culated this is equivalent to the two equations 


! 

Z 

+ R i 
jN+l 

R , + R 

5N-1 5 n 

R i + R i 

iN-2 5N-I 

•• R o 

+ R 1 


u° 

1 

+ q, +1 


0 

R i 

W 

+ R 1 

jN +2 

R l + R l 

iN-2 £N+1 

R 1 + R 1 ’ 
5N+3 jN 

•• R i 

+ R 2 


C ! 

+ c ^ 


0 

r 2 

; Vi 

r i ; Va 

R o + ^-3 

•• Rj 

5N-2 

+ R 1 
5N-I 


C 1 

T ~ 2 

+ C 1 
^N+3 

= 

0 

R 1 

+ I^ 

R o + Vi 

R , + Va ’ 

•• R 

iN-l 

+ R 1 

¥ 


v> 

+ C 1 
^N + 2 


0 

R o 

+ v, 

R , + ^ 

R 2 + Vi • 

•• R 
I N 

+ R ! 

jN+l 


_% 

+ C 1 
^N+l 


I 


(16a) 


and 
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R 0 R N+1 


*N +1 

r n 

r n-i 

R |n+i 

R, 

2 * 


R N ~1 

R N "2 

*• Ri 

R 1 

jN-l 

r n-i 

R N -2 

R N -3 

*• Rj 

? N “ 1 

R 1 

§N “2 

R, 

T v+1 

R l 

Jn 

R 1 

■* R t 

R 0 


“ R 1 
2 N + 1 

R 1 

^N-1 

R ^N + 2 

R 2 

“ r n-i 

- R, 

R 1 

^N-2 

" " 

R 1 

“ R n - 2 

~ R i 

2-N-l 

R 1 

2 N “ 3 

- ft, 
2 N 

‘ R 0 

R N-3 

1 R 1 

R 1 

1 R 2 ■■ 

• Rj 

jN-2 

“ R 1 


R 1 r n 

R 0 “ R N-1 



R 1 r n 

R 2 " R N-1 

R 1 “ R 1 
jN |N+1 


The last equation can be written as 


C 0 + C NU 


I 

C 1 + Cn 


0 

C 2 + C N ^ 


0 

c, + c, 

jN ^N+l 

= 


C 1 

^N+l 



C 1 

jN + 2 



_ C N+i _ 


0 


(16b) 


R, - R, 

JjN £n + 1 

R, - Ri 

J«-| jN+2 

r 2 r n _, 

R, - R. 

2«-l 

Rj - R. 

^N-2 ^N*-l 

Rj ~ R n -2 

R. - R, 

JN-2 ^N-I 

R, - R. 

t N -3 

R 0 - Rn -3 

R 0 - Ri 

Rj - R 2 

R, : R, 

5N-2 1*1-1 


R 1 ~ R 0 " R N + 1 


Cj 


T 



5 N +1 



K o " r n-i r j " r n 


C ^N +2 


0 

r i ” Rn-2 R 2 “ R n -i 


^♦3 


0 

, , - Rj Rj Ri 

*“1 jN jN + | 


C N + 1 


6 


8 » '"V 1 ” 


C 0 + Cn+I 

R " "»-> "■ r |n 


Cj * C„ 

"n-I "n-J • 


C 2 + 

>• i *•_ 


% + *>♦«_ 


where the C 0 + C N+l column matrix would have been found from Equation 16a. Thus, if R n is 
Hermitian symmetric the original system of equations can be solved by solving successively 
two systems each of which is about half the order of the original. 


Naturally, if further conditions are imposed onR n further results may be obtained. The next 
section is devoted to quite special choices for R n with a study of what can be said about C N “ + \ in 
these cases. Before proceeding to these, though, we give one other example. 

Suppose R n = a n k where K is a nonsingular Hermitian symmetric constant matrix and a n is a 
scalar function of n satisfying a_ n = a n in accordance with Equation 1. The fundamental Equation 15 
then becomes 


N + l 


E 

r = n 


a N + l-k-r 


KC 


kO 


I 


which after the substitution X r = KC r is 


for 0 i k % N + 1 , 


N+l 

^ Wk-r x r = s k0 i for o I k g n + i . 

r = 0 

Write X r in the form X r = \ r I +Y r where k r is a scalar function of r so chosen that 

N+l 

°N«-k-r ^ = s «o for 0 * k !■ N + 1 . 

r=0 


(17) 
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The matrix of coefficients of \ r is assumed nonsingular. Then we find 


J2 Wk-r K 1 +Y r) '= 1 YL a N + l-X-r + 

r=0 r=0 r= 0 

N+l 

E 


X N+l-k-r 


which means that 


^ ^kO + ) ^N + l-k-r 


in ^ j. 

E 


^N+l-k-r ^ 


follows that C r 

il 

*1 

1 

Let 


~“n + 1 

“n 

Vi • • • 

a i 

a o 

°K 

“m-i 

“N-2 ' - ’ 

a o 

a l 

“H-l 

“n-2 

v 3 . • / ■ 

a i 

a 2 

a o 

a i 

a 2 

% ‘ 

a N+: 


\+ 


a centrosymmetric persymmetric determinant. Then the system of Equations 17 can be solved for 
\ r , and ^ N+1 has the form 

" + 1 ( X) Vl ' 


Consequently, the mean square error is simply 


C 


-l 
N + l 


(-ir 


EXAMPLES 

Phase Space Process 

Regard x t as a one-dimensional continuous parameter process and denote its derivative by x t , 
The x t process has the Fourier representation (Reference 1) 


/• 00 

Xt = J 

*' -on 


e 2r,ikt dZ(\) 
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so that the formal derivative x t is given by 


/•CO 

x t = 2nik e 2w ^‘ dZ(X) 

J -00 


If we were to sample these at equal intervals we would expect these to go over to 


e 27TiXn dy(\) 


rU 2 

2tt i 

*'- 1/9 


k e 2niXn dy(\) 


Thus we are led to define the two-dimensional "phase space process" X n as 


r x i r 1/2 m— r 1/2 

x e 277i\dy(\) - e dY (M 

L n J *'-1/2 L J *' -1/2 


where 


dY - [^iAj d y( X ) - Ady(X) . 

To obtain the covariance matrix of X n we determine the expected value of dY dY T : 

E{dY dY 7 } - E {a dy A 7 dy } = AA 7 E{dy dy} = P(\)f'(X)S(A.-/u.)d\d/i 

where 


P(A) = AA 7 


X ~2rri\. 
2rriK 47 t 2 A 2 


Here f ' (A) is the spectral density function for the one -dimensional x n process. Consequently, 


r 1/2 

K* m+ „x7} = e^P(\)df(\) 

J -I/O 


e 27 ' imA - dF(\) ; 
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that is, the spectral density function F' (\) of the X m process is equal to that of the x m process 
multiplied by the matrix P(M. It is of interest to determine the counterparts of the above formulae 
in the real case. As seen in the work leading to Equation 5 the one-dimensional process will be 
real if dy = dyj + idy 2 with dy x (k) being an even function of X and dy 2 (k) f an odd. Then, upon 
letting du = 2dy j and dv = - 2dy 2 , we obtain, as in Equation 5, 

J *l/2 

(cos 2jrnk) du(k) + (sin 2jmk) dv(A.) 
o 

and 

r 1/2 

x n = I - 2nk ( sin 2nnk) du + 2nk (cos 2rmk) dv . 


This is representable in the two-dimensional form as 


X 


n 



(cos 2rrnk) dU(\) + ( sin 2rrnk) dV(\) 


(20) 


where 


du 


du 

2 7rk dv 


and 


The du and dv satisfy Equation 7, 


dv 


dv 

2/rk du 


E{du(X) du(/x)} 
E{du(\) dv(yu)} 

but du and dv as given by Equation 21 


- E{dv(X) dv(^)} - g' (X) S(\ -/i) dk dyu , 

3 E{dv(\) du(ytf)} - 0 , 

do not satisfy Equation 7. In fact, 


E{dU(\) dU ( T w) } 
E{dU(A.) dV^j} 


E{dV dV<^)} - Ag' (\) >(A- M )dAd H , 

- E{dV(\) dU TJ ■- Bg 1 (l)>(\-fi)d\d/i , 


(21) 
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where 


A 


1 0 
0 4^ 2 A 2 


and 


' 0 - 7tt\ 

2 TTk 0 


Consequently, 


R E 


r 1/2 /M/2 

{X m+n X n T } = (cos 2n (m + n) k) (cos 2rmfx) E{dU dU T } + (cos 2n (m + n) X ) ( s in 2mi,u) E{dU dV T } 

J o J o 


* (sin 2n (m + n) \ ) (cos 2’'iv ) E{dV dU T } 4 (sin 2 ,r (m + n) \ ) (sin 2”n/i) E(dV dV T } 


f 

J o 


(cos 2?rmK ) Ag' ( A ) dA - (sin 2'rniA ) B g ' ( A ) dA 


or, when written out explicitly, 


r 1/2 

Jo ° 


1/2 /• 1/2 

:os 2"tn\ g' ( v) d\ 2" \ sin 2 T -m\ g' (A) d\ 

Jo 


I" 

J Q 


* 1/2 r \/2 

2n | A sin 2nmk g' (A) dA 4^ 2 I A 2 cos 277mA. g' (A) dA 

Jo 


(22) 


Observe from Equation 22 that one general characteristic of a phase space process is that the off- 
diagonal elements of the covariance matrix are negatives of each other. Also, R m 21 is the derivative 
of R m u with respect to m and R m 22 is the negative of the second derivative of R m n . 

If g(A) is a step function, so that g' (A) has the form 


k 

s' <M = 


j= i 
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then Equation 22 becomes a sum of matrices 


- E’f 


cos 2jTmk- 


2rrk . sin 2rmk ^ 


2nk.^ sin 2miW + 77 2 k ? cos 2rrmk . 


It may be noted that the matrix P(A.) in Equation 19 is singular and so, therefore, is F' (\) . 
Thus, we are not in a position to apply the outlined Fourier approach starting with Equation 12 for 
prediction on a finite past for a phase space process. 

Degenerate Quasi Phase Space Process 

Let the covariance matrix have a value approximately that of Equation 23 with k = 1. That is, 


cr 2 cos ( 2rrvnk ) + ep m 


2nkcr 2 sin 2 mn\ + ep 12 


2nkcr 2 sin 2nm k + ep 21 fyr 2 k 2 a 2 cos 2rrmk + ep_ 


P + ep 

m 1 m 


With p_ m - P„ T , e « 1 . 


We shall show that in this case the mean square prediction error is small; it is, in fact, of the order 


The demonstration makes fundamental use of the relation 


p P _1 P T = p ( 

r n r 0 m r n~m \ 

which follows from a simple matrix calculation. 

We use the equation connecting the prediction error with the covariance matrix in the form 


•^N + l-p^r Cf I § p0 


for 0 $ p i N + 1 


Take the convolution of this with the resolvent kernel 


S „p I " PNM-nP0 _1 Wp 


p - 1 § 
r 0 N+l,p 


for 0 ^ n < N , 0 $ p - N + 1 , 


for n = N + 1, 0 - p - N + 1 , 
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and get 


N+I N+l 


N+l 


y * C np y | ( P N+l-p-r +e ^N + l“P“r) C r ^ ^ C np 1 


P 0 


p=0 r=0 


p=0 


N+l N+l 

y~^ t y^' C np ( P N+l-p-r +€ ^N+l-p-r) 
r= 0 [_p=0 


C r = s n0 I • 


Now, for 0 - n £ N , 


N+l 

E 

p=0 


N+l 


C P 

u np ^N+l-p-r 


y ( S np 1 P N + l~n P 0 1 *N + l, P ) P N + l-p-r 


p=0 

P 


“ P P “ 1 P 

N+l-n-r Jr N + l-n ^ 0 


- p - p p “ 1 p T = o 

r N+l-n-r r N+l-n 0 r u 


from Equation 24. Again, for 0 = n ^ N , 


N + l 

E 

P=0 


C np ^N+l-p-r 


e (^N + l-n-r P N+l-n P 0 1 ^ r ) 


whereas for n = n + 1 we obtain 


N+l 

y \ c n+i,p p n+i 

p=0 


IN f 1 

■ 2 ] p °' 


s p = p _1 p T 

N + l.p ^N + l-p-r *0 r r 


P-0 


and 


N+l 

6 y C N+1.P ^N + l-p-r " eP 0 1 Pt ‘ 
P=0 


Consequently, the result of taking the convolution of Equation 15 with c is 


IN T 1 

0 + € y ^ (^N+l-n-r ” P N+l-n P 0 1 Pr ) “ ^n0 * 

r = 0 

N+l N+l 

2] p 0 “ lp r T C r + P ° _I P * Q * = 0 


for 0 i n N 


for n - N + 1 


r ( 26 ) 
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The first of these equations is of the form ePC = j where P and J are independent of e and the 
second equation upon multiplication by e is of the form e(P + eP) c = 0 with P, P independent of e . 

If e - 0, eC must remain bounded away from zero since the right-hand side of Equation 26 stays 
fixed and nonzero. This implies that = 0(e), if exists (which will be the case for general 
values of p ). 

The fact that the convolution of P N+1 _ p _ r with c np vanishes for 0 ^ n ^ N is connected with the 
fact that the matrix P(A) in Equation 19 is singular. 

Quasi Cosine Process 

Take the covariance matrix to be 


R -* cos 277mA. 0„ + ep 

m ^ U ' m 


with p = p T , e « 1 , 0 < A < 77 , 

' — m ' m 7 


where Q 0 is a nonsingular symmetric constant matrix (i.e., independent of m). At first glance this 
covariance matrix looks simpler than that of the example for the degenerate quasi phase space 
process but the resolvent kernel is more involved. It is, in fact, 


r ( 

|s„p - [sin 2n (N + l-n)\] (esc 2^)8^ + [sin 2n (N - n) \] (esc 2n\) S N+1 p j I 

b np = ^ ( CSC2 S Np “ ( CSC ( COt 2n K ) Q 0 "‘ S N + l,p 


1 ^N+l, F 


for 0 = n ^ N “ 1 , 


for n = N , 
for n = N + 1 . 


The crucial computation is a bit of trigonometry showing that for o = n^N - 1 


N + l 

A 


p- o 


b 


np 


cos 277(N + l“p~r)A 


0 . 


The result of convolving b with Equation 15 is 


r 


iui 

0 + e^~'-|p N+1 _ n _ r - [sin 2n (N + 1 - n) i\.](csc 2 tt\) p y _ t + [sin 2 n (N - n) A.] (esc 2rr\) pf^C T = S n0 I 

0 for 0 = n I N - 1 , 

N+l N+l 

(CSC 2rrK) (sin 27rrA) C r e (esc 2rrk) [(esc 2rrk) p i _ r - (cos 2 t/A) p/JQq” 1 C r - 0 for n = N , 


N+l N+l 

( cos 2rrr\) C r + p r J Q 0 -1 C r 


for n = N + 1 , 
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and reasoning as in the previous example leads us to the same conclusion; i.e., that here also 

c n+! = 0 ( 0 - 


Quasi White Noise Process 

Take the covariance matrix to be 

R m = s m0 Qo + ^ with P _ m = pJ , e « 1 , 

where again Q 0 is a nonsingular symmetric constant matrix. This is a highly degenerate case and 
fairly strong results can be obtained. The equation to be solved is 

N+l 

£ (Wp-r,0 Qo + ^N+l-p-r) Cr = I S P0 for 0 g P * N + 1 . (27) 

r= 0 


If e is set equal to zero it is easily verified that the solution is C r = 0 , 0 i n , C N+1 = Q 0 _1 . Thus, 
if we think of the solution C r to Equation 27 as a function of e , C r - C r (e) , then C r (e) = o(l) , 

0 ^ r ^ N , and C N+1 ( e ) = 0(1) as e 0. Now take the convolution of Equation 27 with the kernel 


5 np 1 6 ^N + l-n Qo 1 S N+l,p 

5 , 


N+l, p 


for 0^n<N, O^p^N + 1 

for n = N + 1 , 0|p^N + l 


and obtain 


[>+ 1 

[^N + l-n-r,0 Qo +£ (^“^r,o) ^N+l-n-r + ® ( £ 2 )] 


8 n 0 1 


for 0 = n £ N 


22 ( S "° Q o + g Pr T K = 0 for n = N+l 


Observe that through the terms in £° and e 1 the first of these equations is of the same form as 
Equation 27 except that the ranges of r and p have both been reduced by one. This is made even 
more apparent by rewriting the equation as 

N+l 

0( £2 ) C 0 + 22 [W* 0 Q 0 + eP N+^-r + °( £2 )] c r = 8 P 1 1 for 1 = P = N + 1 , 

r= 1 
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wherein the transformation n = p + 1 has been made to make the ranges of r and p the same. This 
suggests a reapplication of the above convolution procedure. The convolution kernel of choice now 
is 


^np ^ * € ^N+l-n Qo 1 ^N+l,f 


S N+l,p 1 


for l^n^N, l^p^N + 1, 
for n = N + 1 , 1 = p = N + 1 ■ 


Then upon taking the convolution as before we get 


IN T- I 

O(^ 2 )C 0 + J2 [Wn- r ,oQo +£ (l- 8 :>W-n-r + 0 ( £ 2 )]C r 


S ml 


for 1 I n % N 


IN T* 1 

°( £2 ) C 0 + 2 ^ t S n Q 0 +€ ^l-r +0 ( 62 )] C r = ° fOr n = N + 1 . 

r= 1 


By a little manipulation the first of these can be thrown into the form 


N+l 

O( e 2 )C 0 + 0 (e 2 )Cj + [^ N +3- P -r,o Q 0 + e, °N+ 3 - P - r + 0 ( e2 )J C r = S p2 I for 2 ^ p | N + 1 


We continue in this fashion and ultimately obtain 


rj+i 

0 (£ 2 ) c 0 + 0 (e 2 ) Cj + ■ • • + 0 ( £ 2 )C n + 22 [ S N + N + 2-P-r,o Qo + ^N + N + 2-p-r + °( £2 )] O r = S„ iN+I I for p = N + 1 , 


or 


0(^ 2 ) c 0 + 0(6 a )C 1 + --- +0(€ 2 )C n + [q 0 + 6P 0 +o( £ 2 )]c n+1 = I 
Now as seen earlier C r = o(l) for o^r^N as e - 0. Hence we have 

[Qo + ^o + 0 ( e2 )] C n+1 = I + o(e*) , 
and this may be seen to be equivalent with 


c N ;\ = Q 0 + + o( e2 ) • 


(28) 
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Exponential Decay Process 

Take the covariance matrix to be 


R m = Q 0 - V. . 

where a is positive and real and Q 0 is a real nonsingular symmetric constant matrix. This covar- 
iance matrix is of the form a m K which was investigated earlier in the section entitled ’’Prediction 
Error Matrix Equation. 11 We saw in that case that solving the matrix equation for could be ac- 
complished by evaluating two determinants defined on a m . However, we shall not use that result 
now but instead will obtain by the methods employed in the other examples. The equation to be 
treated is 


F e _ a |N+ i- p -r | Q 0 c r = I s p0 for 0 = p = N + 1 


and the kernel that suffices is 


5 np I - 5 n + l,p e ~ a I 


Vl.p 1 


for 0 = n ^ N 
for n - N + 1, 


In fact, 


N+l 

22 V np e- a l N+1 -P-H 

N+l 

= F ( S np I-S n+ i. p e- a i; 

P~0 

p=0 


- J e -a|N + l-n-r| _ j e ~a-a ]N“ 


= J e ~a|N+l-n”rJ ^ _ e a[N + l-n- 


j j e -a|N + l-n-r| ^ _ gO) 


1 I e -«|N+l-n-r| ^_ e -2aj 


i |N + l— p-r| 


0 = p = N + 1 , 
0 i p ? N + 1 . 


for 0 < n i N 


for N - n - r = 0 , 
for N - n - r < 0 , 


The convolved equation is therefore the diagonal system 

N+l 

22 1 e~ a l N+1-n-r l (1 - e-a«) Q 0 C r = S n0 I 

N + l 

22 1 e "° r Qo C r = 0 


r=N+l-n 


for n - N + 1 . 


0 = n $ N 
0 < n < N 


for 0 = n = N , 
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All we need do is take n = 0, as this equation is merely 


I e° ( 1 - e ~ 2a ) Q 0 C N+1 = I 


which gives 


C N+l = (l-e- 2a )Q 0 . 

Quadratic Exponential Decay Process 

Take the covariance matrix to be 


R = e -an 2 +bn q = n 
m ^0 m ’ 

where a is positive and real, b is pure imaginary, and Q 0 is a nonsingular Hermitian symmetric 
constant matrix. These qualifications are imposed to insure that the covariance matrix fulfills 
Equation 1. This example is most simply treated by applying a succession of convolutions, the 
first being with the kernel 

fs I “ S e - a(2N + l-2n)+b j for 0 < n < N , 

np n+l,p 

U 1 “ “S 

np 

UWp 1 for n - N + 1 . 

This convolution produces the system 


N+l 

L 


r=0 


j e -a(N+l-n-r) 2 +b(N+l-n-r) 


(l - e _2ar ) Q 0 C r 


for 0 = n ^ N , 


N+l 



r = 0 


for n - N + 1 


But the factor l - e 2ar vanishes for r - 0; so we now have a reduced system to deal with: 


N + l 


E 

r= 1 


I e 


-a(N+2~p-r) +b(N + 2~p-r) 


Qo (i 


r )c r 




for 1 5 p S N + 1 . 
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This is of the same general type as the original equation and so we proceed as before, the kernel 
now being 


for l^n^N, 1 ^ p I N + 1 , 

for n = N + 1, 1 = P = N + 1 , 

with the resulting convolved equation 

N+l 

j e -a(N+2-n-r) 2 +b(N+2-n-r) ^_ e “2ar^ ^ C T = 8^1 for 1 = n = N , 

r=l 

N+l 

I e -a(l-r) 2 +b(l-r) ^ _ e ~2ar^ = 0 for n = N + 1 , 

r=l 

which again provides us with a reduced system since 1 - e _2a < r-1) vanishes for r = 1. We continue 
in this way and ultimately obtain 





_ c fi -a(2N + l-2n)+b -r 

° n+l, p C A 


ivri 

^ Ie -a(N+N+2-p-r) 2 +b(N + N + 2-p-r) Q o | 1 „ e -2ar^ 1 _ e -2a(r-l)j ... ^_ e - 2 a(r-N)^ Cr =: J 


for p = N+ 1 


which gives 


N+l 

c N+ i = ri(i- e - 2ap ) q„ 


Note that this is independent of b. 


SUBSAMPLING PREDICTION 

One question of some importance in prediction theory is to what extent the prediction error is 
affected if a subprocess of the given process is used instead of the original process itself. To be 
specific let us compare the linear prediction error for a given multivariate wide-sense stationary 
process x 0 , X, , X 2 , - • • X N+1 = x (K+1)/a with the error for an equally spaced subprocess of this, 
x 0 , X p , x 2(8 , X 3/3 , ■ ■ ■ x (K+1)j8 , the prediction in each case to be of the type X N+1 = X (K+I)/3 on the 
basis of the preceding ones going back to X 0 . According to the linear least squares theory the 
predicted values for X n+1 in the two cases are 

N 

X N+ 1 = - B r X r 

r=0 
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and 


(K+l)/3 


ft 

. ■ -L 


B (K+l)/3 B s/3 X s/3 ’ 


In view of the wide-sense stationarity we have K s/S = R s/3 where R s/S is the covariance matrix of the 
subprocess. The system of equations for C s/5 » B (K+1)yS is then 

R (K + l)yS C 0 + ^K/3^73 + R (K-l)£ C 2/3 + *** + R o^(K+1)/3 “ 1 * 

^K/3 C 0 + R (K-1)^S + R (K-2)/S C 2y3 + **' + R / C (K+l)/3 “ 0 ’ 

R (K-l)yS C 0 + R (K-2)/3 C j 3 + R (K-3)yS^2y0 + ’** + R 2/S C (K+1)^ 0 ’ 


R 0 C 0 + ^J C /3 + R 2£ C 2/3 + *'* + R ( T K+l)yS C (K+l)yS 0 ‘ 

It will be noted that this array is essentially the same as the original array with the exception that 
terms corresponding to the elements cast out in the subsampling process are suppressed. 

It appears that a determination of the difference in the two prediction errors C~ + \ - C^ 1)/8 for 
the general case is not feasible. In special cases certain information may be gained. For example. 
Equation 28 shows that for a quasi white noise process C N ~* - C ( "J 1)/g = o(e 2 ), since the terms in e° 
and • involve only R 0 . Also, the last two examples are simple enough that the error difference 
can be obtained explicitly. In the exponential decay process R s/3 = e” a l s/3 l Q 0 = e~ a ^l s[ Q 0 so all we 

need do is treat the subprocess the same as the original with a replaced by a/3 . Therefore, 

2 2 

-C( K+1 )£ - (e' 2a ^ - e“ 2a ) q 0 , and in the quadratic exponential decay process R s/3 = e~ a P s +b ^ s q 0 , 
whence 


c - 1 - r ~ 1 
Sl+1 (K+l)yQ 


N+l K+l 

TI (l “ e~ 2ap ) - n (l - e- 28 ^ 2r ) 

P=1 v ' r=l 


Qo • 


DIFFERENCING OF MULTIVARIATE PROCESSES 

The first-order divided difference in the multivariate case is defined as 


X[su] = (s -u) _1 {X s -X u } , 
the second-order divided difference is defined by 

X[tsu] = (t -u)' 1 {X[ts] -X[su]} = (t - u)' 1 (t - s) _I (s - u) _1 {(s -u) X t + (u - t) X s + (t - s) X u };(29) 


28 


and so on. Just as in the case of differentiation the differencing of a polynomial gives a polynomial 
of lower degree. In particular, the second-order divided difference of a second-degree polynomial 
(with column matrices for coefficients) is a constant (matrix), and the third-order difference is 
zero. 

Suppose 5 0 , Ej, « 2 , and S 3 are random variables with E 2 and H 3 orthogonal, and let \ be a real 
number. Then the second-order divided differences of the process 


+ 3. t 


H 2 t2 


S 3 e 


27riA.t 


form a wide-sense stationary process; that is, 



+ h, 


s i + h , Uj +h] X [t 2 + h, s 2 + h, 
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is independent of h : 
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and this does not involve h . 

This result extends immediately. The process 




x, = £> f * £>«, e“ v 

k =0 k=l 

for any set of real numbers, \ k , and H 0 , E 1? * * * H p+n random variables with s p , E * * • E p+n ortho- 
gonal has wide-sense stationary p th order divided differences. 

The algebraic polynomial part of x t represents what is commonly called the "trend,” while the 
trigonometric polynomial part is wide-sense stationary. In practice, the trend is commonly re- 
moved by fitting a polynomial to the data, say by least squares, leaving as residue the wide-sense 
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stationary part. As the above work shows a wide-sense stationary process can also be obtained by 
forming the differences of sufficiently high order terms of the original process. Of course, the 
resulting wide-sense stationary process is not the same as the one obtained from removing the 
trend by curve-fitting methods. The former represents the differences of some order of the latter 
(to within the error in the curve fitting). 


(Manuscript received January 21, 1966) 
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